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1.  INTRODUCTION 


Among  linear  multistep  (LMS)  methods,  the  Adams  family  is  found  to 
be  the  most  effective,^  to  this  date,  for  solving  initial  value  prob¬ 
lems  in  ordinary  differential  equations  (ODI-.'s).  Adams  methods  are 
impractical  in  solving  "stiff"  equations  because  prohibitively  small 
step  sizes  are  required  for  accuracy.  To  achieve  equal  effectiveness 
of  Adams  methods  for  solving  stiff  equations,  efforts  were  made  by  Cer- 
taine,-  who  suggested  a  generalization  of  Adams  methods  for  step  number 
k  =  1,2;  by  Bjurel,^  who  modified  Adams  methods;  by  Norsett,'^  who  gave 
an  A-stable  modification  cf  Adams-Bashforth  methods;  and  by  Verwer,^ 
who  generalized  LMS  methods  with  zero-parasitic  and  an  adaptive  princi¬ 
pal  root.  Jain^  developed  A-stable  methods  by  means  of  Uermite  inter¬ 
polation  polynomials,  and  Murphy^  employed  Newton's  divided  difference 
representation  of  the  Uermite  interjiolation  polynomials  to  develop  a 
family  of  A-stable  methods.  Lambert®  introduced  multistep  methods  with 
variable  matrix  coefficients. 

Tliis  report  shows  that  an  appropriate  choice  of  the  characteristic 
polynomial  coefficients  of  the  nonlinear  multistep  (NLMS)  methods  yields 
the  generalized  Adams-Bashforth  and  generalized  Adams-Moulton  methods 
(GAB-CiAM) .  These  methods  are  shown  to  have  portions  in  common  with  the 
methods  of  Certaine,  Bjurel,  Norsett,  Jain,  Murphy,  Verwer,  and  Lambert. 
An  examination  of  how  NLMS  methods  yield  GAB-GAM  follows.  A  section  of 
theory,  describing  the  consistency,  stability,  and  convergence  of  GAB- 
GAM  methods,  will  be  outlined.  In  addition,  the  property  of  A-stability 
of  GAB-GAM  will  be  outlined.  GAB-GAM  can  solve  stiff  equations  effec¬ 
tively.  However,  this  is  not  a  restriction;  GAB-GAM  can  be  used  to 
solve  nonstiff  equations  as  well,  although  this  is  inefficient  for  that 
class  of  problems. 

We  have  included  here  an  application  of  GAB-GAM  to  a  class  of 
underwater  acoustic  parabolic  wave  equations.  A  package  of  ANSI  FOR¬ 
TRAN  (GABM)  programs,  designed  to  implement  GAB-GAM  methods,  is  intro¬ 
duced.  The  important  features  of  the  GABM  package  are  discussed  in  a 
separate  section.  Test  results  of  GABM  are  compared  with  a  set  of 
previously  published  results  by  other  techniques.  Computations  were 
performed  on  a  PDP-11/70  computer,  using  double  precision  arithmetic. 

A  listing  of  ANSI  FORTRAN  (GABM)  programs  is  included  in  the  appendix. 
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2.  FORMULATION 


Initial  value  problems  of  ODl  's,  j;onerally,  can  he  expresseJ  by 


y'  =  yltg)  =  Vg  . 


(2-n 


A  class  of  systems  of  first  oriler  ODF's,  called  stiff,  arises  in  the 
applications  of  chemical  kinetics,  reactor  kinetics,  missile  j5*^idance, 
.  .  etc.,  and  o^ten  takes  the  form 

y'  =  Ay  ♦  g(t,y):  y(tg)  =  yg  , 


where  Re{X(A)}  <  0  and  eigenvalues  of  A,  differ  greatly  in 
magnitude.  We  choose  to  deal  with  equation  (2-2),  since  the  applica¬ 
tions  frequently  appear  in  this  form. 

In  solving  equation  (2-1),  k-step  LMS  methods^  can  be  applied  that 
take  the  form 


i  =  0 


a  •  N' 

»■  n  +  i 


C-.^) 


where  0,  |ag|  ♦  |Bgl  >0. 

LMS  can  solve  equation  (~-l)  effectively  when  ||:>f/3yl|,  tlu'  norm 
of  the  Jacobian  matrix,  is  small.  Our  approach  is  to  generalize  equa¬ 
tion  (2-3)  such  that  when  ((9f/3y(|  is  large  LMS  can  be  used  to  solve 
equation  (2-2)  effectively.  11\is  led  to  what  we  refer  to  as  NI.MS  meth- 
ods^O  that  possess  the  form 


Ah(k-i) 
a .  e  '■  y 
1  '  n>i 


'I’l  ■  (Ah)g 
k  1  ^  > 


i=0 


(2-4) 


where  ^  0,  jag)  ♦  h (-l^i^g (Ah) )  |  >  0. 

In  carrying  out  the  complete  development,  A  is  assumed  to  be  a  con 
stant  nonsingular  matrix.  Methods  to  handle  A(t)  will  be  explained  in 
a  later  section.  It  can  be  seen  easily  that  when  A  =  0,  equation  (2-2) 
reduces  to  equation  (2-1)  and  equation  (2-4)  reduces  to  equation  (2-3) 

*0  '^ki^^^^  ”  *^i'’  identity  matrix.  This 


where 
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shows  that  NLNiS  methods  include  LMS  as  u  subset.  Tlie  theory  of  the 
NLM.S  methods  also  generalizes  the  theory  of  the  LMS  methods,  as  has 
been  demonstrated  in  reference  10. 

The  complete  formulation  of  NLMS  involves  the  accurate  determi¬ 
nation  of  i>j^j(Ah)  with  the  selection  of  a;.  Tlxis  appears  in  its 
entirety  in  references  10  and  11.  A  brief  outline  of  NiXS  methods  is 
given  below. 

We  first  define  the  NLMS  operator. 


My(t):h|  =  ^  ♦  i! 


-  h  ^  4'l^i(Ah)g(t  ♦  ih,  y{t  ♦  ih))  . 


Then,  for  a  constant  matrix  A,  we  express  equation  (2-2)  by 


(e'^V) 


e  g(t,y),  yCtfl)  =  yq 


(2-b) 


and  integrate  equation  (2-6)  over  the  interval  obtain 


...  r  ^n*i 

^  *J 


e''^^n^i‘*'^g(t*,y)dt'  . 


(2-7) 


If  we  expand  y(tn  ♦  ih)  and  g(t',y)  of  equation  (2-7)  in  a  Taylor 
expansion  around  t^^,  substitute  into  equation  (2-5),  and  simplify,  we 
obtain 

ik  j  00 

^  ^  C^(Ah)g^j^  ,  (2-8) 
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where 


C^(Ah)  = 


If'"!/ 


(t'  -  t^)-'dt’ 


■  Z'T!"' 


In  view  of  the  first  condition  of  consistency, 


i  “i  =  f  • 


which  implies  that  {  }  of  equation  (2-8)  is  equal  to  zero,  we  obtain 


Llytt);hl  =  ^  C^(Ah)g^-’^  . 
j=0 


(2-10) 


Since  A  is  assumed  to  be  a  constant  nonsingular  matrix,  we  can  mul¬ 
tiply  both  sides  of  equation  (2-9)  by  A-’^^,  and  we  obtain 


pAnetn^i-t  _  tn)Jdt’  = 


UAIO" 


(2-11) 


Using  equation  (2-11)  and  simplifying  equation  (2-9),  we  obtain 


i=0  £=0 


(2-12) 


We  require  that  Cj(Ah)  =  0  for  j  =  0,1,2 . .  but  C  i(Ah)  ^  0. 

This  requirement  defines  an  NLMS  method  of  order  p.  The  usual 
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consistency  conditions  for  LMS  are  general ized . This  permits  formu¬ 
lating  NLMS  methods  in  a  matrix  form,  leading  to  the  solution  of  <f>i<i(Ah3 
explicitly  when  the  characteristic  polynomial  coefficients  are  selected. 
In  requiring  that  Cj (Ah)  =  0  for  j  =  0,l,2,...,p,  we  also  select  the 
order  p  to  coincide  with  the  step  number  k;  this  selection  ensures  a 
unique  solution  so  that  ♦j^jCAh)  can  be  solved  explicitly. 


We  now  select  the  characteristic  polynomial  coefficients  =  1, 
=  -1,  0^.2  =  =  •••  =  oq  “  This  selection  leads  to  the 

formulation  of  Adams  and  generalized  Adams  methods. 

For  a  constant  nonsingular  A,  we  multiply  equation  (2-12)  by 
(Ah) let  1|a|]  -*-0  for  j  =  0,1,2,. ..,p-l,  and  define  (0)  = 
as 


S  “i  ■  JT  ^  i’H  .  (2-13) 


i=0 


i=0 


Equation  (2-13)  relates  to  oj  and  gives  the  second  condition  of  con¬ 
sistency  of  LMS  methods. 

When  i^k)((Ah)  =  0,  equation  (2-12)  is  expressed  as  a  predictor,  as 
follows: 


I 


I  +  (k  -  l)Ah 


I 


I  +  kAh 


V  -  DAh]"'  V  (kAh)'" 
^  m!  ^  m! 

m=0  m=0 


-e 

I 


Ah 


(2-14) 
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’Ah 

0! 

O  ■ 

1  1  •  •  •  1 

♦  kO  (''^*') 

0  1  •  •  •  (k  -  1)1 

'i'kl  (Ah) 

• 

• 

• 

• 

o 

(Ah)l' 

(p  -  1)! 

• 

0  1  ...  (k  -  1)P'^I 

• 

♦k.k-l^Ah) 

Kc  can  find  ♦j^j(Ah)  as  follows: 


Substituting  <Ji|,-(Ah)  of  equation  (.2-15)  into  equation  (2-14),  we 
have  the  Generalized  Adams-Bashforth  (GAB)  method. 


When  (♦i.i.CAh)  /  0,  equation  (2-12)  is  expressed  as  a  corrector,  as 
follows: 


I 


I 


1  ♦  (k  -  l)Ah 


I  ♦  kAh 


V  (kAh)™ 
2.  ~in! 
m=0 


I 


(2-16) 
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1  -(Ah)'' [-(Ah  -  -  I],  -(Ah)'‘[-e  ♦  (I  ♦  Ah)] 


-(Ah)'^[-(1  -  ♦  (I  ♦  ^)], 

-(Ah)'^[(21  -  (Ah)-)e^^  -  2(1  ♦  Ah)]. 


-(Ah)'^[-(I  ♦  ♦  (I  ♦  |Ah  ♦  (Ah)')]  . 

-(Ah)'‘*[(I  -  i(Ah)2)e^‘'  -  (I  ^  Ah  .  j(Ah)')]. 

-(Ah)  ^[-(31  ♦  Ah  -  (Ah)')e  +  (31  +  4Ah  +  4(Ah)')], 

-(Ah)''’[(3I  ♦  2Ah  -  4(Ah)'  -  (Ah)^)e^^  -  (31  +  SAh  +  3(Ah)')], 
-(Ah)''*[-(1  +  Ah  +  i(Ah)')e^‘^  (I  +  2Ah  +  ~(Ah)'  +  (Ah)^)]  . 


Substituting  ♦|^j(Ah)  of  equation  (2-17)  into  equation  (2-2),  we 
have  the  Generalized  AJaros- Moult on  (GAM)  method.  Simplifying  the  terms 
of  ♦i  -  (Ah)  of  GAB-GAM  and  letting  |  1a|  [-*■  0,  one  obtains  the  Adams-Bash- 
fortn  and  Adams-Moulton  methods,  respectively. 
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3.  nnioRY 


Details  of  the  comi'lete  theory  of  NLMS  methods  with  regard  to  con¬ 
sistency,  stability,  and  convergence  can  be  found  in  reference  1('.  The 
theory  of  NLMS  methods  is  automatically  applicable  to  CAR-CI.-XM  methods. 

A  sketch  of  the  theory  of  (iAB-tiAM  methods  Is  given  liere. 


CONSISTENCY 


An  NLMS  method  is  defined  to  be  of  order  p  when  Cj (Ah)  is  selected 
to  be  0  for  j  =0,1,  .  .  .,p,  but  Cp^2(Ah)  ^  0.  This  condition, 

Cj (Ah)  =  0,  implies  that 


1  im 
h->0 


k 


MS 


Ah(k-i) 

.  e  ^  \ 

1  '  n  +  i 


k 

i  =  0 


0  . 


(3-1) 


This  equation  yields  the  consistency  for  CAB-GAM  methods. 

If  we  let  il-All-*  0,  we  have  shown  that  equation  (2-12)  reduces  to 


s 


TJ 


k 


i  =  0 


(3-2) 


This  is  an  alternate  confirmation  of  the  second  condition  of  consis¬ 
tency  for  the  LMS  methods. 


STABILITY 

The  characteristic  polynomial  of  NLMS  methods  is  defined  by 

p(X,c)  =  e^''^p(0  ,  (3-3) 

where 

k 

P(C)  =  ^  a.c'  .  (3-4) 

i  =  0 
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In  view  of  the  choice  of  a.,  GAB-6AM  are  strongly  stable.  This 
choice  also  gives  the  first  condition  of  consistency  of  both  LMS  and 
NLMS  methods. 


CONVERGENCE 


Tlieorem 

GAB-GAM  methods  are  convergent . 

This  theorem  has  been  proved  in  reference  10. 


A-Stability 


Dahlquist^^  defined  a  method  to  be  A-stable  if  the  numerical  solu¬ 
tion  IITj^II  0  asymptotically  as  n  ->■  «•  for  the  differential  equation 
y'  =  Ay,  where  Re{X(A)}  <  0.  GAB-GAM  methods  have  this  property. 


Theorem 


GAB-GAM  methods  are  A-stable  in  the  sense  of  Dahlquist. 

If  we  apply  GAB-GAM  methods  to  solve  the  problem  y'  =  Ay,  which 
implies  g(t,y)  =  0,  the  GAB-GAM  methods  produce  the  exact  solution  to 
the  problem;  i .e . , 


At 


nAh 

=  e  y. 


lim 

n-K» 


We  calculate  the  e*^^  by  Fade  approximation. 
11/^11  establishing  the  A-stability, 


If  Re{A (A) }  s  0,  then 
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■>.  imim.i;mi:ntati()n  oi-  cab-ham  Mirnioos 

'»h™~ 

axe.  namoly.  variablo-step-si ’e'  .>rol  .f  V"  *  "^'‘’n’orat  od  in  the  pack 
starting  teclniii,ue.s.  and  the  hind  in^! 

A(t).  The  ijsaee  of  fhi«  ^  ^'”''ahle  and  t  i me -dependent 


A(t).  The  usaee  of  t  is  J  “''X  variable  an, 


i me -dependent 
A  descrijition  of 


PK0CR/\M  STRUCTURi; 

diajjram  below.  *  mon}>  the  programs  is  shown  by  the 


START 


INVhRT 


Ibe  main  control  proi-ram  (MAIN)  processes  ill  m,  > 
a  secondary  control  subroutine  (Oll-m)  iVnm  . 
procedure;  calls  for  the  CAB-CAM  netln  i  - 

■step-size  changes,  predict^  i.ul  sol  f-starter . 

snb.iect  to  the  user-rei,u  ired  tol er  ince  'u  convergence, 

tor  printout.  The  self-starter  ISTMrn'  r\  I  ^  Prepares  the  results 
matrix  inversion  (  INVPRT)  md  vit  li  reJi’  '‘I’proximat  ion  (PADi:). 

t.ne  formats.  CPN  must  Ik^  supph^^i'by^■IH\,.^;";V‘'' 

l■■|;ATllRl;.s 

la  ollaiiKial  hy  lialviny  or  Jouhliiio  lo,',  *■  I'ooiloil.  rlio  stop  s‘i;o 

-tor  I,,:  is  „so.r  t„^„lio  t  f;  ‘ 

. . . . '  -r;p:  r'::rys^n^;:*s;::; 
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2.  Prc\l ict-aiul -Cor root -iii-Timos,  This  procoduro  is  omployod 

when  a  var ial)le-step-si2c  tecTuiiquc  is  used.  Setting  m  =  .>  allows  the 
corrections  to  be  carried  out  uji  to  three  times. 

‘'^*-^li~!^tart  ing  Procedure.  I'irst  order  generalised  /Vlams-Bash- 
forth  method  is  used  as  a  self-starter.  Tlu'  step  size  used  for  the 
starter  is  10  times  smaller  than  tltc  user-suggested  steji  size. 

4.  A  is  a  I'unction  of  t.^^  In  solving  stiff  equations  of  the  form 

y'  =  A(t)y  +  g(t,y)  ,  (4-1) 

the  progr;un  is  constructed  to  decompose  the  equation  into 

y'  =  A(t.)y  ♦  (Alt)  -  Alt^ly  +  g(t,y)  .  (4-2) 

so  that  the  N1>1S  methods  c<n  be  applied.  Since  the  above  equation  is 
stiff,  RelA(A(t))}  <  0  for  all  t.  If  |  |A(t)  -  A{tp||  can  be  main¬ 
tained  .small  enough,  the  above  problem  can  be  solved  accurately,  ilio 
decomposition  is  ]ierformcd  automatically,  provided  the  indicator  lAT  is 
set  to  1 .  As  one  can  see,  after  the  decomposition,  the  new  g(t,y) 
becomes  {A(t)  -  A(t^)}y  +  g(t,y);  only  the  original  g(t ,y)  needs  to  be 
defined  in  tl)e  CRN  subroutine.  The  UMAX  is  recommended  sucli  that 
I  |A(t)  -  A(t.)||  <  hP.  An  examiilc  is  given  in  the  section  of  numerical 
examples . 


INPUT  INFORMATION 

A  listing  of  inputs  is  given  below  with  tlieir  definitions: 

ERR  User -required  tolerance 

T(l)  Initial  time 

TMAX  Final  time 

IPC  =  0,  fixed-step  size 

^  0,  variable-step  size 
INDEX  =  0,  explicit  methotls 
/  0,  implicit  methods 

KSTEP  Step  number  (also  order  of  the  method  <  4) 
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thea 


H  Initial  step  size 

(MAX  Largest  step  size  allowed 

HMIN  Saallest  step  size  allowed 

N  Order  of  the  system 

YZERO  A  vector  array  to  store  the  initial  values. 

Soae  default  values  arc  built  in,  in  case  the  user  failed  to  supply 


They 

are 

UMAX  = 

0 

.1 

D- 

■0 

imiN  ° 

0 

.1 

D- 

-5 

ERR  = 

0 

.1 

D- 

■7 

H 

0 

.1 

D- 

1. 
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S.  iV  'A'  -  ■‘.V  ADF.ftWJTF.R  VaVR  ’RCPAGATICN 

Th^  ic'^isf-c  v^ve  ^r.'caear.on  "Totsieins  lar.  :e  jxrr“ssec 

1  i  1  y  -)y  jri  ell.o'^AC  ^nuation, 

A  O  *  =  ;  ,  5-L 


p  th<^  ico'AStiC  ores-^ure, 
<  i  <;  ?h<^  v^v^Tiiimh^r ,  anrl 


n  n^r,z,  rhe  .nd«x  of  refraction. 


•»s<;'jmi fiij  the  i nhomoi^snei  t . e*;  are  alowl/  yar^-'in?  in  ran^e,  a  cyL^.o- 
'Ificall/  a/n>wetrjc  ^eometr/,  and  a  hartnonic  source,  eauaticn  3-1  car. 

transformed  into  t>>o  uncoupled  parabolic  vave  eouations.  Ine  para¬ 
bolic  equation  stands  for  the  transmitted  field,  representing  the  eave 
propagating  in  the  range  direction.  The  other  parabolic  eouatccn  stands 
for  the  reflected  field,  representing  the  .iiave  propagatcng  cn  the  depth 
direction.  Because  of  the  slowl/  varying  property  of  cnhcmcgenectces 
in  range,  the  reflected  field  is  negligible.^  Therefore,  one  needs  only 
to  sol'/e  the  transmitted  field,  and  Tappert^  ^  did  this. 


Tappert's  algorithm  is  Vnown  as  the  split-step  Fourier  algorithm 
and  is  considered  to  he  the  only  effective  technique  for  solving  the 
parabolic  wave  equations  in  the  underwater  acoustic  field,  with  some 
limitations.  Por  one  class  of  problems,  GAB-CAM  methods  can  produce 
reasonably  accurate  solutions  rapidly  without  restrictions. 


first,  let  us  take  the  ODK  approach  to  treat  the  solution  of  para¬ 
bolic  wave  ef(uations.  We  obtain  a  .sy.stem  of  ODE's  where  the  second  par¬ 
tial  derivative  with  respect  to  the  .space  variable  is  discretized  by 
means  of  a  second  order  central  difference,  as  follows. 


A  parabolic  wave  ef|uation  in  one-space  dimension  takes  the  form 


II  =  a(kn,r,z)ii  ♦  b(kn,r,z)u 


zz 


(5-2) 


Ilf  scret  i  z i ng  iij,^  by  a  second  order  central  difference  brings  equation 
fri-2|  to 


dll  b 

m  m 


r  -  =  fl  ti  *  -  2u  +  u  ,1 

ilr  mm  n?*  m+1  m  m-1^ 


(5-3) 
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Equation  (S-3)  is  a  system  of  first  order  ODE's  that  can  be  put  in 
the  matrix  form 


Jui 

3r 

2b, 

“  ’IP- 

bi 

0  •  •  • 

0 

0 

Ul 

bl 

dU; 

IT 

*>2 

2b, 

'iiT- 

^2 

-4  *  ■  * 

0 

0 

U2 

0 

• 

dr 

■ 

0 

0 

0  ... 

m- 1 

b  1 

m-1 

• 

♦ 

• 

• 

“m-1 

h2 

'^m- 1 

0 

- 1 

i 

0 

0 

0  ... 

2b 

m 

“m  -  h2  _ 

“m 

b 

m 

h2  m*l 

•  > 

where  Uq,  u  are  two  boundary  points  at  range  r.  To  put  the  system 
above  in  a  short  matrix  form,  we  get 

u'  =  A(r,z)u  +  g(r,z,u)  .  (5-5) 

In  selecting  a  range  step  size  for  solving  equation  (5-4),  using 
GAB-GAM  methods  with  a  constant  matrix  A,  one  will  choose 

h|lUlllUkk(Ah)l|  <  1  (5-6) 

for  corrector's  convergence. 

In  the  range -independent  case,  g  can  be  decomposed  again  such  that 
||3g/3u||  =  0.  This  implies  that  a  large  step  size  can  be  used.^l 
Note,  also,  that  the  ODE  approach  formulates  the  problem  without  other 
restrictions. 

A  first  order  GAB-GAM  is  used  with  the  predictor-and-corrector 
feature  to  solve  the  following  problem.  The  equation 


ui 


-  l)u 


♦ 


(5-7) 


describes  a  shallow-water  wave  propagation  in  the  range  direction  in  the 
region  indicated  in  figure  1. 


14 


TR  601 


Figure  1 .  Region  of  Wave  Propagation 

The  source  is  placed  at  9.8  m  below  the  surface  and  the  bottom 
depth  is  19.5  m.  The  bottom  is  rigid;  thus,  a  Neumann  boundary  condi¬ 
tion  is  assumed.  We  examine  the  wave  propagation  up  to  400  m.  The 
sound  speed  profile  is  described  in  figure  2. 


Figure  2.  Sound  Speed  Profile 

The  initial  sound  speed  Cq  =1.5  km/s,  frequency  =  25  Hz,  the  sur¬ 
face  condition  is  u  =  0,  and  n(r,2)  =  c(r,z)/cQ. 

The  solution  was  obtained  through  four  different  methods:  normal 
mode,!'*  explicit  finite-difference,  first  order  Adams  PC^,  and  the  first 
order  GAB-GAM.  The  normal  mode  solution  is  very  suitable  for  this  prob¬ 
lem;  therefore,  we  use  it  as  our  reference  solution  for  comparison  of 
accuracy.  Solutions  by  explicit  finite-difference,  Adams  PC^,  and  GAB- 
GAM  all  agree  almost  exactly  with  the  normal  mode  solution;  the  differ¬ 
ence  is  in  the  speed,  as  tabulated  below: 
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h 

Time 

Explicit  finite-difference 

0.35  X  10-5 

8  m  1 1  s 

Adams  rC^ 

0.34  X  10-** 

10  m  9  s 

GAB -GAM 

1280 

23  s 

The  step  size  is  that  which  is  required  to  give  an  accuracy  of  10"'’. 
What  makes  CAB-GAM  so  efficient  is  that  the  resulting  OPIi  of  equation 
(5-7)  takes  the  form  u’  =  A(:)u  +  K(r,z.u).  In  this  application,  g  is  a 
function  of  u  and  ||3g/3u||  =  |  |h|^/(Az)^  |  |  =  [  [ Co/2iif (Az) -  |  |  ,  which  is 
not  a  small  number  for  reasonably  small  (Az)  .  However,  the  rigid  bottom 
allows  us  to  restate  the  problem  in  a  desirable  way.  Notice  that  the 
last  equation  of  the  system  is 


(Az 


-yU  ,  ♦  (a 
)  2  m-1  \ n 


\  b 

m  )  m 

(Az)-A  m  (Az)-“ m+1 


0  . 


(5-8) 


The  rigid  bottom  condition  enables  u^^j  to  be  expressed  as  u^; 
therefore,  equation  (5-8)  can  be  rewritten  as 


(Az) 


m-1 


(  .  AJ  - . 

V"'  (Az)v  m 


(5-9) 


This  reformulation  defines  g  to  be  0;  therefore,  ||3g/3u||  =  0;  this 
means  that  a  larger  step  size  can  be  used.  In  fact,  the  solution  to 
this  problem  by  the  GAB  method  has  a  simple  explicit  form. 


•'^n+l 


Ah 

=  e  y 

^  r 


(5-10) 


Hence,  subject  to  the  accurate  computation  of  e^^,  the  solution  can 
be  found  in  one  step.  We  can  see  that  for  a  range- independent  matrix  A 
and  rigid  plane  parallel  boundaries,  GAB-GAM  methods  have  definite 
advantages . 
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6.  NUMERICAL  EXAMPLES 


GAB-GAM  methods  are  designed  to  solve  stiff  ODE's  in  the  form 
y'  =  Ay  p(t),  where  p(t)  is  a  low  order  polynomial  in  t  and  A  is  a 
constant  matrix,  in  the  absence  of  roundoff  and  Fade  approximation 
errors.  In  the  following  examples,  we  present  a  semistiff  system,  a 
nonlinear  stiff  system,  and  a  stiff  system  with  A  =  A(t).  The  last  two 
examples  are  not  in  the  format  recognizable  by  NLMS  methods.  We  will 
show  how  these  problems  can  be  brought  into  equivalent  forms  such  that 
GAB-GAM  is  applicable.  The  starter  employed  is  a  tirst  order  GAB  with 
a  step  size  10  times  smaller  than  the  initial  step  size  user  suggested, 
and  all  the  examples  are  started  at  time  *  0.  For  every  fixed-step 
size,  only  one  matrix  inversion  is  needed  for  the  algorithm.  An  addi¬ 
tional  matrix  inversion  is  needed  for  the  Fade  approximation.  There¬ 
fore,  only  two  matrix  inversions  are  needed  to  solve  the  whole  problem 
for  a  fixed-step  size. 


EXAMPLE  1,  SEMISTIFF  SYSTEM^O 
Problem 


Exact  Solution 


Eigenvalues  of  A 


{1.  -10}  . 


Results 


Method 

(Order) 

*max 

Step 

Size 

Number 

of 

Steps 

Number 

of 

g  Calls 

Number 

of 

Inversions 

Error 

GAB (3) 

10 

10 

1 

3 

2 

41  x  10-»2 

GAB 

0 

Exact 
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Remarks 


NLMS  methods  are  designed  to  solve  this  type  of  problem  exactly, 
subject  to  the  accurate  computation  of  eAh.  in  this  case,  it  is 
expected  that  NLMS  methods  will  produce  accurate  results  with  large  step 
size  without  requiring  large  computational  expense. 


EXAMPLE  2,  NONLINEAR  STIFF  SYSTEM^S , lb . 17 
Problem 

y[  =  -0.013y2  -  1000y^y2  -  2500y^y3  ;  yj(0)  =  0  , 
yl  =  -0.013y2  '  1000/1/2  J  /2C0)  =  1  . 

y^  =  -2500yjy3  ;  yjCO)  =  1  . 

Reference  Solution 


yi(48)  =  -0.194533895680D-5  , 
y2(48)  =  0.611474831446  , 
y3(48)  =  1.388950571516  . 

The  above  problem  is  in  the  form  y'  =  f(t,y);  yCtp)  =  y^.  This 
problem  is  not  in  the  form  directly  acceptable  by  NLMS  methods.  Since 
the  system  is  autonomous,  we  can  easily  restate  it  appropriately. 
Transform  the  dependent  variable  such  that  Zj(t)  =  yj(t);  Z2(t)  =  y2(i) 

-  1;  ZjCt)  =  y3(t)  -  1.  Consequently,  z(0)  =  0.  We  now  expand  f(z) 
about  the  critical  point;  i.e., 

z’  =  fCz)  =  f(0)  +  zf'iO)  +  O(z^)  =  f'(0)z  +  g(z)  , 

where  A  =  f ' (0)  =  (Sfi/Bzjjj-o  and  g(z)  =  f(0)  +  0(z2).  The  matrix 
thus  found  is  singular.  The  value  of  a33  is  adjusted  to  -0.01  to  keep  A 
nonsingular  while  similarly  readjusting  g(z)  .  Therefore,  we  obtain  an 
equivalent  system,  as  follows. 
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■-3500 

-0.013 

o" 

’^r 

-0.013  -  lOOOz^z^  -  2500Zjz’ 

^2 

= 

-1000 

-0.013 

0 

"2 

-0.013  -  1000ZjZ2 

.^3. 

-2500 

0 

-0.01 

/3. 

-2500ZjZ3  +  O.OlZj 

with  Zj(0)  =  Z2(0)  =  ZjCO)  =  0. 

Eigenvalues  of  A 

1-0.00993,  -0.01,  -3500.01307}  . 


Results 


Method 

(Order) 

t 

max 

Step 

Size 

Number 

of 

Steps 

Number 

of 

g  Calls 

Number  of 
Matrix 
Inversions 

Error 

GAB (2) 

48 

16 

2 

A 

9 

00 

X 

lo-** 

GAB(2) 

X 

o 

10-** 

Jeltsch 

Remarks 


For  autonomous  systems,  an  expansion  about  the  critical  point  of 
the  right-hand-side  function  produces  the  desired  format  for  GAB-GAM 
methods . 


EXAMPLE  3,  STIFF  SYSTEM  WITO  A  »  A(t)5.8,18 
Problem 


-60  +  0.125t 

10 

r0.125t1 

0 

y'  = 

0.2 

-0.2_ 

♦ 

1 - 

o 

1 _ 

;  y(0)  = 

_o_ 

Reference  Solution 

yj(lO)  =  2.344886  x  10-2  ^ 
y2(10)  *  1.301528  X  10-2  , 

yiC400j  =  27.110701  , 
y2(400)  -  22.242211  . 
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Introducing  a  new  variable,  >’3  =  t,  brings  the  above  problem  into 
an  autonomous  system. 


0 

1 

1 _ 

10 

0.125' 

0. 125yjy 3 

0 

>■’  = 

0.2 

1 

0 

0 

y  ♦ 

0 

:  y(0)  = 

0 

0 

0 

-1 

1 - 

♦ 

1 _ _ 

0 

■■  * 

Eigenvalues  of  A 

{-0.2335,  -1,  -59.9065}  . 


Results 


Method 

(Order) 

^max 

f 

Step 

Size 

Number 

of 

Steps 

Number 

of 

g  Calls 

Number  of 
Matrix 
Inversions 

Errc 

ir 

10 

1 

53 

68 

X 

10-11 

GAM (3) 

GAM(3) 

7 

49 

X 

10-^ 

Lambert 

400 

1 

397 

2783 

50 

X 

10-3 

GAM (3) 

44 

X 

lo-** 

Lambert 

Remarks 


NLMS  methods  are  developed  for  constant  A.  For  A  =  A(t),  the  prob- 
lejfl  is  restated  in  a  form  acceptable  to  NLMS.  The  main  strategy  is  to 
construct  g(t,y)  to  be  a  low  order  polynomial  in  t  or  a  slowly  varying 
function  in  t.  Note  that  the  solution  at  t  =  10  is  almost  two  orders  of 
magnitude  better  than  Lambert's,  despite  the  large  step  size.  The  com¬ 
parison  at  t  =  400  given  by  Lambert  may  itself  be  in  error  since  he  used 
Runge-Kutta  methods  to  achieve  yC400). 
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7.  OTHFR  ADAMS- LIKE  DF.VrLOPMRNTS 


CliRTAINE'S  DF.VFL0PMI;NT“ 

Certaine  considered  the  differential  equation  in  the  form 

>■'  =  -uy  +  g(t,y)  ,  (7-1) 

where  D  is  a  diagonal  matrix.  Multiplying  both  sides  of  equation  (7-1) 
by  e*^^  and  integrating  it  over  the  interval  [t^.t^],  Certaine  obtained 

'■a  , 

e^f'''-^2)g(f,y(t’))dt’  .  (7-2) 

1 

where  h  =  t2  -  tj;  g(t',yft'))  is  approximated  by  a  polynomial  in  t', 
call  it  a(t'),  of  degree  n  +  1;  and  n  is  restricted  to  be  ^  2. 

The  integral  of  equation  (7-2)  can  be  written  as 


r^2  ,  " 

/  e^^''''’'2ig(t'.y(t'))dt'  ^  2  ' 

■^t,  j  =  0 


where 


A  (t  -  ^l)-’ 


g(t,y(t))  »  a(t)  =  y  -  -  —  -j 

^  j  !  li 
j=0 


and 


C.  =/"  'eD(t'-t2) 

3  /  j!h3 


It  is  immediately  seen  that 


C.  =  D-l(I  - 


(7-3) 
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and 


C 


n+1 


For  the  first  order  method. 


Xj  =  ♦  (Cp  -  C,)gj  ♦  Cjg^  .  (7-4) 

For  the  second  order  method, 

>2  ■  c-'"V,  .  (c,  -  ^,)s„  .  (C„  -  2C;)8,  .  (^j  •  C,)ej  .  (7-5) 

If  we  take  A  =  -D,  simplify  (Cg  -  Cj),  and  substitute  it  into  equa¬ 
tion  (7-4),  we  get 

♦  h{-(Ah)-2[[(I  -  Ah)e'^*’  -  I]gj  [(I  +  Ah)  -  ] } .  (7-6) 

which  is  exactly  rhe  first  order  (.'..AM  method.  Further,  it  can  be  veri¬ 
fied  that  Certaine's  method  of  order  2  is  the  second  order  GAM  method. 
Therefore,  we  see  that  the  GAB-GAM  methods  include  Certaine's  work. 


BJURtL'S  MODIFICATION^ 

Bjurel  considered  the  differential  equation  of  the  form 

Pj^(d/dt)y  =  f(l,y)  ,  (7-7) 


where  P  is  a  polynomial  of  degree  n  with  constant  coefficients.  Bjurel 
generalized  LMS  methods  and  obtained 


0;  a^(h)  =  1  . 


i  =  0 


where  a^(h),  6^(h)  are  mesh  dependent. 
Bjurel  defined 


k 

P(i;.h)  =  ^a.(h)!;'  . 

i=0 


(7-8) 


22 


TR  6011 


k 

o(c,h)  =  ^  e.Ch)?^  . 
i=0 

In  order  to  give  rise  to  the  Adams-like  methods,  Bjurel  first  consid¬ 
ered  pCc.h)  is  of  the  form 


V  Xih 

pCC.h)  =  n  ^  ,  (7-9) 

3  =  1 

where  Xj  are  roots  of  Pj^  =  0,  and  then  selected  the  remaining  (k  -  n) 

roots  of  P„  =  0. 

N 

The  resemblance  to  the  ordinary  Adams  formulas  can  be  seen  by  put¬ 
ting  N  =  1  and  Xj  =  0. 

NORSETT’S  method'* 

Norsett's  A-stable  modification  of  the  Adams-Bashforth  methods  con¬ 
sidered  the  problem 


y'  =  f(t,y);  y(to)  =  y^  (7-l0) 

over  the  interval  I  =  [tjj,b]. 

Choose  t^el  and  define  t  .  -  t  =  ih.  Write  y'  =  f(t,y)  as 

y’  +  Py  =  Py  +  f(t,y)  .  (7-11) 

where  P  is  a  constant  matrix.  Approximate  T(y)  =  y'  +  Py  with  a 
Lagrange  interpolation  polynomial  at  points  i  =  0,1,..., q. 

Integrating  T(t)  between  t^^  and  t^^^  and  using  methods  of  gener¬ 
ating  functions,  we  obtain  the  following  formula: 

q 

>■„.!  ■  ‘  "  I  *  "''n-l)  • 

i=0 


23 


whero 


m /"  l’l>s/-s\  ,  t  -  t,. 


q  =  0.1.:....;  i  =  0.1.....1,. 

It  is  oasily  seon  that  if  wo  take 
and  1’  *  A.  wo  j;or 


0  =  0.  i  =  0,  i;  =  f  *  py 

- 1  II  - 1  •'  n  - 1  ’ 


>n.l  =  ‘'''’’Vn  *  ''(AlO-'n  -  o^'’lf„  . 


17-1.-S) 


This  is  idontical  to  tho  first  order  CAK  methods. 


JAIN'.S  Ml  niou.s^’ 


Usinj.  a  similar  starting;  point  as  Norsott.  .Jain  developed  A-stahle 
iithods  tor  .Stitt  GUI  's  based  on  llermite  interpolation  polynomials,  l.et 


(t  -  a,)lt  -  aj 


'  ‘  -  «j.iHt  -  a.^j)-  .  .  .  (t  -  a  ) 


e.(t)  =  - - . _ L_  j  -I''  "j  +  r  • 

J  Uj  -  aj)(a,  -  a,).  •  •  (a.  -  lr7[)(a. 

be  the  Lagrange  interpolation  polynomials. 


lyt)  =  [1  -  (t  -  a.)e}(a.)]t‘(t) 


(7-15) 


(t)  =  (t  -  a. )e^ (t) 
.)  .1  .1 


be  Hermitc  interpolation  formulas.  Write  y'  =  f(t.y)  as 

y'  +  Py  Py  +  f(t.y)  . 


(7-l(>) 


(7-17) 


Approximating  >■'  ♦  Py  hy  llermite  interpolation  polynomials  yields 
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y'(t)  +  Py(t)  =  ♦  Pyp  >  ^  ’’^i^ 


1  =  1 


i=l 


where 


J _  ,.(2r) 

'n  (2n)! 


K..  =  F^^‘^(C)w2(t)  . 


P(t)  =  f{t)  +  Py(t)  . 


and 


iT(t)  =  (t  -  Bj)  (t  -  a,) 


•  •  • 


(t  -  a  ) 
^  n 


Integrating  equation  (7-18)  from  t^  to  t^^^j  gives 


-Ph 

Vl  =  ^n  *  ^ 


-Pt 


n+1 


S  (II. F.  tl.F!) 

^  \  I  11 


i  =  l 


*  ^  * 


where 


..  i 
.,■/ 


n+1 


n+1 


Pt 

e  h.(t)dt  . 


(t)dt  , 


n  (2n) ! 


/■ 


e  r  pt  fonl  o 

•  e'^F^^"J(C)Tr2(t)dt  , 


P  is  an  approximation  to  -(3f/3y)^  (which  is  chosen  to  be 


P  =  - 


^n  - 

^n  -  ^n-l 


for  a  scalar  equation), 


(7-18) 


(7-19) 


(7-20) 


(7-21) 


(7-22) 


(7-23) 


(7-24) 
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'’ii 


ri/.  1  i 

•  •  •  ^n-l' 


\ 


y'  -  y‘ 

'n  ^n-l 


(7-25) 


for  a  system  of  equations,  and  P  is  considered  to  be  a  diagonal  matrix. 

Equation  (7-20)  can  be  rewritten,  by  changing  variables  and  letting 
s  ”  (t  -  t^)/h,  as 


where 


''n.i  ■ 


f  (H^F,  .  HjF!) 


Lia  1 


♦  R_ 


(7-26) 


„  I  Phs.  , 
Hi  ”  /  e  hj(s 

J r\ 


(s)ds  , 


(7-27) 


rr  {  I’l'Sir 

H.  -  I  e  h. 

'  J  ‘ 

^  t\ 


(s)ds  , 


(7-28) 


m 


^^/'n2(s)ds 

■'o 


♦  0(h^"^^)  .  (7-29) 


Therefore,  the  method  (7-20)  is  of  order  2n. 

The  integrations  involved  to  determine  Hi,Hi  arc  of  the  form 


.1  .  N 


•'o  'i-1  / 


ds  , 


(7-30) 


N  »  0,1, ... ,n. 


ii 
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Take  a  simple  case,  n  =  1,  and  we  have  hj^(t)  =  1;  Fj(t)  =  t  - 
iT(t)  »  t  -  hj(s)  =  1;  Fj(s)  =  s;  and  ii(s)  =  s.  Then 


and 


1 

^  p2h2  ’ 


(s)ds 


1^ 

6  • 


Select  the  coefficients  of  Hj  to  be  0,  and  equation  (7-26)  becomes 

♦  he’’’^|(Ph)-l(e’’’’  -  1)|f^  .  (7-31) 

Define  -P  =  A  and  and  equation  (7-31)  is  identical  to  the 

first  order  GAB  methods. 


MURPHY'S  DEVELOPMENT^ 

Murphy  proceeded  again  from  the  same  starting  point  as  Norsett  and 
Jain.  Murphy  wrote  y'  =  f(t,y)  as 

y'  >  Ly  =  Ly  ♦  f(t.y)  .  (7-32) 

Multiplying  both  sides  of  equation  (7-32)  by  e^^  and  integrating  from 
'n  Vl 

r^n+1 

y^Vl^  *  ^  e^^f(t,y)dt  .  (7-33) 

^n 

Murphy  employed  Newton's  divided  difference  representation  of  the  Her- 
mite  interpolation  polynomial  of  degree  2r  to  treat  the  integral. 
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A  change  of  variable,  letting  s  =  (t  -  t  )/h  in  the  integral  of 
equation  (7-33),  gives  " 


n+1 


e^^f(t,y)dt  =  he'^*' 


/ 


n+1 


e*"^*V(t^  +  sh,y(t^  +  sh))dt.  (7-34) 


The  function  f(tj^  ♦  sh,y(t^  +  sh))  is  approximated  by  Newton's  divided 

difference  formula,  using  the  values  t  ,t  ,t  ,,t  ..•••,t  ,,t 

^  IS  n’  n’  n-1’  n-1  ’  n-r+1’  n-r+l 

Therefore, 


f(t  ♦  sh,y(t  ♦  sh))  =  f  ♦  shf  +  s^h^f  , 

'  n  ^  u  ’  ’  n  n,n  n,n,n-l 

♦  (s  ♦  l)s^h^f  ,  ,  + 

n,n,n-l,n-l 


(7-35) 


.  (s  .  l)2sV*-'f  ,  R  . 

n,n, • • • ,n-r+l ,n-r+l  n 


The  error  term 

R  =  ls(s  +  l)***(s  ♦  r  -  ll^h2nf[t 


n’  n 


n-r+1, n-r+1 


The  f^  in  equation  (7-35)  is  the  i-th  Newton's  divided  dif¬ 

ference  given  fiy  the  quotient 


f  •  -  f 

e  _  1  +  1,  •  •  •  ,1+j  i,  •  •  ‘.i  +  j-l 


and  f  „  is  defined  to  be  equal  to  f ' (t  ). 
n,n  ^  '■  n-^ 

The  Murphy  method  takes  the  form 

y  =  e  ^^y  +hraf  +a,f  +a-f  ,+••• 
n+1  ■'n  ^  0  n  1  n,n  2  n,n,n-l 

^  &  f  1 

2r-l  n,n,n-l ,n-l ,  •  •  • ,n-r+l ,n-r+l '  ’ 


(7-36) 


where 


wi  -Lh  {  _ 
i;  =  h  e  IP. 

t\ 


,  ,  Lsh  , 
(s)e  ds  , 
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Lh 


>  (-1)^  — p  -  y  (-1)^  -i — , 

.pi  jTo  (1-'')^ 


(7-37) 


and  the  P^"^(s)  can  be  calculated  by  the  recursive  formulas  below: 

=  ft  *  m  -  l^D('^) 


2m 


m  =  1,2, •••,r  n  =  0, 1 , 2, • • • , 2m-l; 


,(n) 


l,(s)  =  (s  >  m  -  • 


p-^(s)  =  (s .  m  -  i)p^;;!i(s)  *  np^;^:;^s) ,  (7-39) 


(n-1), 


(7-38) 


m  =  1,2, •••,r,  n  =  0, 1, 2,  •  •  • , 2m;  with  P^°^(s)  =  P  (s)  =  1,  p'-'^^fs)  =  0. 

0  0  m 

The  choice  of  a.  =  0  for  j  ^  1  produces  equation  (7-40)  from  equa¬ 
tion  (7-36),  which  coincides  with  the  first  order  GAB  methods  for  j  =  0. 
This  choice  gives 


.(-o, 


y  ,  =  e 
n+1 


a,  =  (Lh)-l(l  -  e'^N  , 

+  hj(Lh)-l(l  -  e'^‘')]f^  . 


0 

■Lh 


(7-40) 


VERWER'S  FORMULATION^ 

Consider  an  LMS  formula  of  the  type 


^"♦1  ~  ^^t^^n'^n^^n-ft-l  *  *^n  ^  ®t^*^n'^n^^n-^l-t  ’  (7-41) 


8,=  1 


)l=l 


where  is  a  matrix  that  equals  or  approximates  J(t^,y^),  the  Jacobian 

matrix  3f/3y.  Based  on  approximate  choice  of  rational  function,  Verwer 
arrived  at  a  multistep  formula  with  zero-parasitic  roots  and  an  adaptive 
principal  root  from  equation  (7-41).  This  formula  leads  to 

k 

Pn.l  =  P(Vn»'r,  *  "n  2 ’  Vn.l-tl  <’-«) 

t=l 


Equation  (7-41)  can  be  put  in  the  form 
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where 


^n+1  *  **^n*l-*,  *  ' 

fc»l 


(7-43) 


p  =  y |A„(h  J  )  +  h  B„(h  J  )M  ,  „}  ; 
^  I  n  n  n  n-*  n+l-lj 


!l=l 


If  we  use  a  Pade  approximation  for  P  and  let 
equation  (7-43)  is  identical  to  the  GAB  method.  Note  that  Verwer  does 
not  consider  implicit  techniques  at  all. 


LAMBERT'S  METHODS® 


To  solve  the  problem  y'  =  f(t,y);  yCtg)  =  yp,  consider  a  class  of 
multistep  methods  with  variable  matrix  coefficients,  as  defined  below: 

k  k 

•  "“itVlVi  ■  "Ilei  •  •  <’•«> 

i=0  i=0 


where  y.  =  y(t.);  f.  =  f(t.,y.);  t.  =  t^  ♦  ih;  |a.(t)|  <_  a  constant;  and 
|bj(t)l  constant  for  every  ^clt^.t^^^]. 


A  class  of  the  above  methods  has  been  considered,  as  follows: 


I 


s 

k 

s 

y  .  =  hT 
•'n+i  ^ 

S=1 

i=0 

S  =  1 

(7-45) 


s 

and  llQpll  f.  q  and  l|aj^®^l  +  nonsingular. 

s=l 


To  label  a  class  of  the  above  methods  as  Method  1,  Lambert  selected 

1;  k  -  1;  p  -  1,2;  aj°^  =  1;  aj^^  =  -1;  bj®^  =  1/2;  bj®^  =  1/2; 
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^  =  0;  and  a^^^  =  0.  This  selection  gives  the  first  order  Adams- 
Moulton  method.  A  class  of  these  methods  is  used  to  solve  problems  in 
examples  2  and  3  of  section  6  with 


Q 


n 


3f 

3y 


=  t 


n+2 
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8.  CONCLUSIONS 


NLMS  methods  have  the  advantage  of  allowing  the  use  of  a  large  step 
size  for  solving  stiff  equations;  CiAB-GAM  methods  have  the  advantage 
that  they  are  strongly  stable  and  require  the  least  number  of  operations 
in  the  class  of  NLMS  methods. 

A  good  implementation  of  CAB-CAM  methods,  certainly,  will  solve 
initial  value  problems  of  stiff  equations  effectively.  CAB-GAM  methods 
are  not  restricted  to  solving  stiff  equations  only,  as  example  1  of  sec¬ 
tion  6  shows.  However,  the  CAB -CAM  methods  are  not  generally  competi¬ 
tive  with  ordinary  Adams  methods  for  such  problems. 

Pade  approximation  plays  an  important  role  in  the  approximation  of 
e^*' .  A  number  of  formulas'^  have  been  given  that  approximate  eAli  effi¬ 
ciently.  The  P.AOE  subroutine  listed  in  this  report  is  for  eigenvalues 
being  negatively  large  and  differing  greatly  in  magnitude.  For  mild  or 
semistiff  problems,  other  Pade  formulas  may  be  substituted. 

CAB-GAM  methods  are  designed  to  handle  OOP's  in  the  form  y'  =  Ay  + 
g(t,y),  where  A  is  a  constant  and  nonsingular  matrix.  For  A  =  A(t),  a 
decomposition  can  be  made  to  transform  the  problem  into  a  form  accept¬ 
able  to  GAB-GAM  methods. 
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Appendix 


A  LISTING  OF  COMPUTFR  PROGRAMS 


A  listing  of  ANSI  FORTRAN  computer  programs  follows.  Among  the 
programs,  the  matrix  inversion  subroutine  is  not  listed.  The  user  can 
supply  a  matrix  inversion  from  his  computer  library.  The  main  program 
describes  the  solution  setup  of  example  2,  from  section  6  of  the  main 
text.  This  problem  can  be  used  as  a  test  case. 


th:s  Pi.-.PM-  conr.jst?  «  St'  Of  s.i*.  tj  :*»TcsMTt 

•  t'-sTf  ctf  o»  The  tcpn 

rv(  T  l.-OT'AV  ♦  0«T,V  I 

Owl«  T-t  CLCStO  ;"«TEB\.#>E  (Till. »"•».. 

OPTIONS  I  t.M'l.tllE -STEP- Size  iPIoEB  STf*  Sltt  At 
SttP-ST«*TINO.  P»t3i:T0B-C0««CT0«. 

*  1$  •  rtTK'ION  Cf  T 

INTuT  PPOANCTCPS  MNVC  TMC  roiLOWINC  OtFiNiTlONtI 

*  •  a-OINENSIONPL  Mltav  PA'fltx 

CO*  U$f*-SPPLtCP  COOO*  TOICRPNCC 
9CrMfl.r  tPLUE  •  *.19-1 

N  INITIAL  STCf  $m.  tCFAULT  VALUC**.**! 

»«IAII  THC  AAXIRU*  STCF  SIZt  TWC  UtC*  ALLOWt  AttfA* 
TO  COIFSIOC*.  KfAUlT  UALUC*I.M 

lAT  tCT't.  IF  A  IS  A  constant 

SCT'I.  IF  A  IS  A  ruNCTION  OF  T 

lom  SCT't.  IF  OtT.V)** 

SCT'l.  IF  ecT.V)  not  •  0 

INOO  SIT«0.  ASCINO  FO*  CMAllCtr  fVTMOM 
SCT'l.  ASCINO  FO*  INALICIT  NCTNOM 

INC  KT'O  FO*  FNIOICTO*  0*  CO««Ccro* 

SCT*t  FO*  FIICCICTO«-ANB-CO*«Ccra* 

*mr  tre*  nurk*  (Lett  than  ai 

A  NUNK*  OF  COUATtONt 

T  l-tlNCNSION  AMAV.  Till  CONTAINS  INITlAl  TINi 

tram  final  TIAK 


•  VZC*0  l-BlfCNSIONAl  AAOAV.  CONTAININC  INITIAL  WALUKt 

tttaitttttiiitiiiitiKttttttiattitaitttMittMttattmtMaMttta* 

CONNON  Ailt.ltl.TIA) 

BinCNSlON  V(4.I*).VZC*0(|*),VNCIItt0l 
^aic  mClSION  A.Cm.N.H^.T.TNAK.V.VlSIlO.S.NniN 
DOuaiC  FACCISION  X.Rl.M.Xl.zi.b.Ml.ui.M 
BOOtLC  FAECISION  NUtXS.VNCH 
SATA  CR*/.tO>S>'.NAAX/l.M/ 

Zl  ••«.  I  B4S3It9SMtD>S 
2C«  «.SUa4T4B3l44SM 
23*  I.SUtSOSTlStSDO 
N*3 
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KSTC^-t 

tPC** 

INOCH't 

ierN«t 

nM)l«4t.Df 


9t 

«t  t  .3  I*#.  £4 

*(2.t  <*-ie««.c* 

•>a,a'»-*.*i3o« 

•  13.1  >*-35M.D* 

•  <3,2»*.M 

•  <3.3)*-*.*lB* 

!•  M>4.C*tH 

T< I  )««.0« 

VZCRO)  1  (•t.M 
VZC40(a)«».M 
vZCROOt't.M 
iriH  .CT.  sm 

C«Ll  DlFEOtCm.N.HMAX.KSTCP.N.TIMX.V.vltRO.ICFN.IPC. 
■  iNOCX.IAT.Z.VnCW.MUSCO.HniNI 

wirccs.tt  MuscD.z 

I  ro«n*T(i*x. -M  •‘.cts.a.iM.'T  •'.cis.s//> 

I  F0«tWT(IX.3(D24.1C.2Xt> 

3  roanMTC/iSK.'MM  ML  CWO«  •'.M4.t»/) 

X2>yt«U(3)*I.D* 

WRITKS.a)  VMCUU  ).Xl.l(i 
URITC(S.2)  Zl.Zi.Z3 
Wt  •(•tS( I VMCW( I )-Zl IrZl > 

W2«e«BS((Mt-Z2)/Zai 

u2*0KAXi(yi.iizt 

y3*0«ISl(K2-Z3»'Z3) 

W3>0nAXi(y2.iai 
WIITC(S.3>  W3 
00  TO  10 
IW 


S  TPO  TiNt  5Ta3Ti«,<.<.  M.H,  y.yN.VOLD.  !▼.  t>»T  1 

ccni^.'H  «( !•,  n  « 

DiniN5:o><  vt 4.ioi.vN« l•t,m.D(4,l•l 

couru  eNicisioN  •.m.t.v.vm.vold 

tAT«  l*,ll.O.|/ 

00  I  ('{.N 
t  VOLDll.l  ••Vll.l  I 
DO  10  JM.KSTCA 

CALL  NLASIll.M.mO.H.VN.lO.II.IT.tAT.IOl 
DO  &  IM.N 
voLD<t.n*v»»<n 
S  V(J«1,|I<VN<1) 

to  COATtNUC 

KTURN 

(ND 


SUtOOt/TINC  AFNTIA.N.T) 
OIPtNSlOA  A<t0.tt> 
DOUDU  AOCCISlOA  A.T 
OCTUON 
CMD 
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SUl'IOtriNE  WMIO.S.H.y. J.T.(*) 

DOUtlE  PNCCSStOn  A(t».t«I.V(4.1*l,Cn«>.T<4l,H,TM 
6i2)«-«.«l3M-|tM.MtYU.t  >«V(J.2) 

J,| )»V(J.3> 

C(l  *  02) 

CO)'.  SD-tlV<  J.  J)  ♦  Th 

RCTUON 

CND 


su»«io<niNE  1)W4T«(*.h,«ns  ) 


I  rWTSI*  INVEBSICD  SUPOUllt*.  CftLlEO  IV  P«UE  OP  HIPS 
«  A  CCPTAJnS  'ME  OPJOIMAL  ELE^'-PTS  AND  REPAIRS  UMALTEPEO 
•  APS  COWTAIPS  tme  Ati(-l) 

I  this  SET-UP  IS  USING  UMOA-:  IIAI  PaThPR  CAi$Tl)«G  DOUlU 
I  precision  6AUSS-J0RDAN  REDUCTION 

I  This  PROCRAfl  IS  REPLACtAlLE  lv  THE  USIP 


POUILE  PRECISION  A(  t*.  !•  I .  AN3(  > t .  !•  ) .0(2  I 
DinENilON  JCdP) 

DATA  "(R/1*/,NC/1R/ 

OID'i.D* 

DO  t  I'l.N 
DO  I  J'l.N 
1  ANSd.JI'Ad.J) 

CALL  DCJRlANS.MR.NC.N.N.PDEX. JC.O) 

IFinOEX  .EO.  1)  GO  TO  II 
RETURN 

1*  URITE  (4.21 

a  FORPArON.aaHHRTIlIM  IHuERSION  ERROR) 
RETURN 
END 


S  PP'  'INf  Pm'-Eim.-.,)' .  B.C  .N  1 


cal:  .Li.’c 
CAL.T  P. 


•  fvrot."NT;«L 

N."')  SU8^CU'l'«l 


By  fade  afppoxINATION 


N'l 

lr(x‘,o)?T»N  -  .101)  3. 21, 21 


D'iceLE  fpt  cii  :0N  P<l«.l«'.P<lP.;|l.i)l|.I»).C<ll,  )•  ),COL'  II I 
P<.>’.,Pif  FfcfoiS’ON  m,ha^E,  v»H)RA 
HA..  I  -M 
00  ?  Id.N 

DO  t  J'l.N 
l(l,J)>I.OI 
ClI.Jl'I.DI 
rd ,  j  )'|.D« 

1  CONTINIjE 
COM  f  )-I.DA 

2  continue 

DO  1?  I'l.N 
O')  16  J'l.N 

cot  (I  t'CCL*  I  )«l>AlStA(J.  I  )) 
l«  CCNTINilf 
17  continue 

XNCPr.COM  1  ) 

DO  18  I-I.N 

irotH^  .CT.  CCLO)  GO  TO  II 
RNoTr-cnLc I > 

It  CONTiml, 

*  roLJ-N  fioei  IS  USED  TO  SEE  UHEThER  t*F<Ai  NEEDS  REDUCTION 


■  ExPlA)'d-I.S«A*Alt2d2.llll(-l  )ld*C.SlR*«ll2/ia.|l 

3  DO  I  I'l.N 

DO  S  J'l.N 

DO  4  R'l.N 

Pd.  J  i'P(I.Jl*Ad.KIIA(R,JI 

4  CONT  tHuE 
Cll.Jl'tPd.J  'tHdZ.MI-Ad.Jl/B.MlIN 

5  CONTINUE  n: 

cd.i  )'fd.i>»i.ip# 


.vy 


A 


>  A 


.NX 


AAA 
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tmS  PAGE  IS  BIST  3’OALITY  lHACTIOAiUtl 
TOOll  aofi  fU«UlSM*D  TO  DUC  . 
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cofWN 

COUlLE  PRCCISIOM  A.AMOim.CM.C.H.MHM.HniH.T.TCA.TIMK 
OCtlU  RStClSIOW  TZERO.V.VN.VNCU,VQLO.V2t*O.MU$C» 

DATA  iZERO.IOAC/t.l/ 

LATO 
JA-t 
MUSCO'H 
DO  4t  1>|,M 
4*  Vl|,l )«VZCR0(tl 
ITtU-O 
ISTEA-D 
m40*TH  ) 

40  ITTR.ITtl»»l 
IH'A 
lAIM-O 
M  CCr«TIf*OC 

(  C'.-AUJ<<TE  A(T*  AT  T-TiO*  ir  lAT  NOT  *  • 

*  IF  KSTCP  CT  I,  CALLS  START  B 

IF(IAT.ME.B>  CALI  Af NTIA.N.TM  M 
01  IFUSTCP.EO.l  .and.  index. CB.Ot  GO  TO  M 
CALL  STA«T(rSTEA,M.H.V.yA,VOLD.IG.IAT) 

M  TCA*T(n«tBLE(FLOA'(«STEA»»tH 

IF(TEA.GT.TnAX)  M«  TAAX'T(n)/DBLE<FlOAT(KSTCA>) 

IFITEA.CT.rnAXI  GO  TO  01 

lf(IM  .CE.  32767)  IM-a 
DO  61  J*1,KSTCP 
ft  T(j4t)*T(J)»M 
lA«».KSTEA*i 
DO  62  JM.IAA 
DO  62  IM.N 
02  VOLDtJ.I  )<v(J,t ) 


B  riXED-STEP-SIZE  IWO,  JMDE<-B  PAECICTOR 

B  IU»6.  INOEX'l  CORRECTOR 

B  OARIABLE-STEP-SIZE  lU  NC  B.  PREDICTOA-CORRECTOA 


IF(10.E9.6.AND.lNOCX.E3.t>  CALI  HUASUSnP.n,  rCLD.N,  VN,  IZE»0, 
6  IH, IG.IAT.JH) 

IftIV.CO.O.AnD. INDEX. HE. t>  CALL  NLnS<XSrEP,H,vOLD.N,VN. lONC. 

6  IN. IG.IAT.JH) 

tPdV.NE.f)  CALL  NLROUSTCP.H.VOLD.N.VN.lZCRO.IONC.ia.IAT.JH) 
Of  00  Of  {■t.N 


bS  'OLPi  CTEs*!  ,  I  '.'.n:  I  ■ 
;f. IU.se. 0  CO  TO  65 
irrLNT  .C).  i  >  CO  TO  69 
o:  *>3  s  •;  .N 

68  VNt  u.  I  '•  VS(  ;  I 
GO  T.)  S2 


t  CORRECTOR  P’  .NCS*  CORRECT  3  "INES  « 

t  S‘’E“-S;3£  CMrtSGEO  BV  A  FACTOR  C^  2  B 


CO 

7#  ICCRPMCCRRTl 

IFiICOKR  .LE.  LNT)  CO  TO  7$ 
H'M.a.SB 


JN-# 

iF<H  .IT.  MAIN)  N'HNIN 
IFIM  .LT.  MAIN)  lNtN-lNIN*l 
IFCniN  .CT.  3l  WRITE  (4.117*) 

117*  FORNATOX.JSnh  reached  HNlN.  NO  CONwCRGENCC  POOSIBLC) 
IFdNlN  .CT.  3)  STOP 
T(1 »*rzE»o 

TEA«T( I  ) 

DO  71  IM.N 
?1  VII.II'VZEROfl) 

CO  TO  SB 

70  CALL  NLRSirSTEP.H.VOLD.N.VNCW.IONE.ICORR.IC.IAT.jH) 
IFILNT  .NC.  I)  CO  TO  70 
CO  TO  82 
70  COnTTTRjE 

ANCRn-*.S* 


88  58  8 
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C 

c 

c 


c 

c 

c 


atiiitittfittKiititttiiitiittittttiiitiiiititsaMMtifitatiitittt 
a  Tc$r  co^fttcTon  a  coNvcauNcc.  usinc  hm  nomi  a 

aataiiatitiiatitiiaiitaaatataaaaaaaaaaaaaaaaaaaaaaaaaaauaaaaaaaaa 

•0  M  J»I,H 

C«.'  ■•(VOLC<KSTCP*l.J)-VICU«J)  I 

.6T.  DalS(«(JI>l  80  TO  M 

•K«n.»«as<a«  j)  I 

cchtihuc 

ir(M«r>ai-cmi  la.u.tiu 

DO  ai  j*t.N 

VCLO((STEP»t,J)*VTCuiJ» 

60  TO  7a 
00  83  IM.N 
V(KSTt**t.l  I.VNCUdi 

atiiitttiKittitiaataaaaaaiiataitaaaaaiiaaaaaaaaaaaaaaaaaaaaaiaaaa 

!i)!Jl’JJ!****************”*””*****»”*”»*»«»«**««««MDM*»«M 

Am'wnaTLN 

NUbCO-OnNaKN.NUSCO) 

DO  as  j't.Ksrcp 
DO  8$  1*1. N 
V(J.I  I » 

iF(j.ca.t)  vztmii fvij.! I 

continuc 

JN«| 

T(t  KTca) 

TZCaO'TCtl 

iF'tu.Eo.a  .NNo.  iNKK.to.ii  irr>i 
IFdU.EO.a  .AMD.  INKM.CO.I)  INDC<*8 
IFCLHT  .18.  II  lM>t 

iFiOADaircA-rfiMi  .it.  <  .id-th  nctudn 
iPdW  .to.  8)  00  TO  88 
Td  1*780 


t  1  •  !  N 

vrfi-o  1 1*  ’  «f*ip*i, : » 

88  ''  '  I .  I  !•>  >KStEP*l .  I  > 

C  tt««tttiMt«it<ti«(««tiaitt(t*«tta«*tittttitiaitt(t««t(atttiittata 

c  t  N'5:H''mN  s'jc;iS4Fn.  m  cohstaniiv  roa  j  tines  befooe  douDiino  a 

C  a  I •••*•**•• ttttitattatttatttaftattataatatttiitaiaaaitttaaattaaaasa 

IS’t«>dSTEP*l 

iri ISTCP  .LT.  31  80  TO  40 

ISTtP.O 

H*d.O<»aN 

triH  .8T.  MANX  I  N«NMa 
ir«0ABS«Tt*-T1IMI  .U.  I.ID-TII  OCTUON 
80  TO  40 

cno 

8 


SUDSCl'TI'-f  M' r?  .“.  V  .  N,  vN ,  t ‘'rt*,  IS,  IT  lA*  Jm  ' 

atttttitiiittitittiiittiiiiiitititiittitiaiiitMtiiiKtiatiia**** 

i  NONtltN.aK  <1  irisTEP  A'.iCP.’ *X-’i  i»a.  G**  I 

»  calle:  dv  eirco  .a  st^at 

1  AaOiaitNTS  ALSEAOV  CEFIAFO  IN  Na;n  PasCPAH 
a  this  paocPAP  calls  3  Si-Iaoutines 

a  INViaTIAH.N.Pl > 

a  6FN1C.....A) 

a  PADC(A, . . ..N) 

a  SOLUTION  LCCTCa  is  stoned  !N  VNdi 


CONRON  AdO. 181.7)41 

5!S'**!2?  *'<«l#.l8).A*«d8.l8).AN3d8.10l,AM4)l8.l8).ATd0,ID) 

•!52fIS2  ^•^••‘•’•^•'••>*»>#).0Hi)4,  ja.iai.uNiTdO.ioi 

DOUtLC  PWtCISICN  FC,P1,AI,A2.A3 

8  JH  AND  IS  APt  INOtCATOaS  *  **  *** 

8  POOCaAA  DOES  INITlALUAd)|N  WMCN  IS*t  AND  JN*0 

8  raosaAN  calculates  and  saucs  an.  CXPIANI.  a  INWKOOC.DMI  fn 

8  IS. 87. 1  ADOVt  CALCUIATIONO  AMt  l«OA$KO 
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lUt  u  NOAtllK  r’lUUliAi^ 
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i 

f 

! 


i 


|F«JH  .CT.  •!  CO  TO  KtTt^ 

PC  2  I'l.N 

DO  2  J'l.N 
A2(I,JI>«.D0 
Md.JfO.CO 
«M(  l,J)>Nt«(I,J) 

M42ii.j»«a.Df 

a  M44(I,J>>O.CO 

ir«N-n  T.i.T 
■  Aid, I  ■•I.Ct/AMII.I  i 
A2d,ll«*ld,tltAl<t.t) 

AHed.l  )-«H«l,|  I 

A3d.t  •■ABd.lllAKl.t ) 

AH3d,|  I’MUC  1,1  IBAH(|,1  t 
CO  TO  9 

7  CALL  INVCirT<AH,N.At» 

DO  21  I'^N 
DO  21  JM,N 
DO  21  KM,N 

AH2(l,J)*A»42(I,J)*AMII,KlfAM(K,J) 

11  A2II,JI«A3M,J)*AI(I,C>IAI<K,JI 

iriKsrcA  .CO.  i >  co  to  f 

DO  22  I'WN 
DO  22  J>t,N 
DO  22  KM,H 

AH](  I ,  J  )>AN3d  .  J  )«AN2(  I,«  ItAMIK,  J  > 

12  A3d,Jt«A3(l,J>*A2d.KI(At(K,J) 

IFIKiTCA  .CO.  21  CO  TO  D 

DO  2t  I<1,N 
DO  23  Jd.N 
DO  23  Kd.N 

13  AM4II,JI«AH4(l,JI*A3d.K)«AIIII.JI 
•  IF(N-il  ID. It. ID 

tl  CAMd.tl'OCKDIAHtl.tn 

SO  t:  14 

ID  CAL'.  PaO£iA.h,FAH,C2.**«,13Ah,C,HI 

14  IF  IS  .ST.  I  1  CO  TO  (IDD,2*D,3M1,  KSTCA 
DO  I  IM.A 

DO  6  J-l.N 
FCU.I,J)*D.PD 
FCt2.l,.’>>D.rt 
FCi3,l,J)«D.DD 
Ald,J)«D.DD 
UMTtUJt'.D.DD 
C2AM(l,JfD.DD 
C3AHtl,J)«D.DD 
I  CONTI hit 

I  UNlTd.t  |.|.M 
DO  3  Id. 4 
DO  3  J*I,N 
DO  3  K«I,N 
FI(I,J.KI*D.DD 
AMI<|,J,KI*D.DD 
3  DNI«  l,J,K  fD.DD 

60  TO  dDD,2»D,3DD)  .  CSTCA 
IDD  CONTInlC 

Ittf  ttlttltttlltlXIttlttllllBtllBttftasll 

t  NCNLINCAR  AULTISTCA  STADTS  NCDC. 

•  tCCINNING  SCCTION  DOCS  INITIALIZATION 


DO  132  Id.N 
VN( 1  ••D.tD 
AmI<2,|,|  mD.OD 
132  AN|i3,I.I  CD. CD 

IFtlS.CT.I  .AND.  INDCX.CD.Dl  CO  TO  131 
DO  ID3  I'UN 
DO  ID3  J>1,N 

113  PltI,JI>-CANI|,J)«UN|Td,J) 

•  1ST  CRICA  CAD 

•  00  icrp  Its  calculates  rNi(i,D> 

I  LOCP  109  C*  IID  CONPUTCS  FINAL  VIN41 ( 

IFlINCex  .NC.  Dl  GO  TO  163 
ID4  |F< IT  .CO.  Dl  GO  TO  IDD 
DO  IDD  IM.N 
DO  IDS  J*t.N 
DO  IDD  K*t.N 
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IM  COMTINUC 

t3t  CALL  CrN(Q,M.N.V,KSTlP.T.AI 
DO  IM  IM.M 
•0  IM  J«I,N 

IN  VNI 1  >•  VN(  I  l*CIINI  I.  J  )«V(  I.  J  >*«MrC« I .  I,  J  )M( i i 
«0  TO  ITS 
IM  M  II*  {•I.N 
DO  II*  J>I.N 

11*  VN(II>VN(I)«CM«(1.J»VI1,J> 

00  TO  ITS 


a  1ST  OVDC*  OAA 

a  conpurc  PHtii.*).  PNiii.i ) 

a  PINAL  aCSULT*  AMt  CALCUlATCD  IN  LOON  lU  0«  IM 


IM  IPIIT  .CO.  *>  00  TO  107 
M  IM  1*1. M 


DO  161  J*l,6 

OJ  I'M 

1*2  r  1 1  1  .  I..’  MC  It  I ,  I  .  J  I  ,(  xlMMtk.,  J  ) 

^m,i.j)-'’i(i.i,j>-Aui,  J) 
lei  Pl«2,I.J>»'>llI.J>»H»All.J» 

16*  CONTImC 
17*  K2>LSrCP«l 

IPlIT  .CO.  •>  02  TO  107 
DO  169  I'l.N 

DO  164  K*l,Ka 
CALL  C^NtO.H.H.V.r.T.Al 
DO  16S  J*1.N 

IPU.CO.l )  NhIO.I.I  )•PHI< 9,1.1  )*P III. I.J  ltd  Jl 
IM  IPlK.Ca.2)  AMI(2,t,|)*ANlta,I,|)*Fl(i.l.JitdJI 

1*4  NNI(9.l.ll*AHU].|.l»«NMl(a.l,it 

169  CONTIN't 

DO  166  I'l.N 
DO  186  K'l.N 

ICC  VN(I  l•VN(l  )-MtA2(|.K  ItPMIO.K.l  >*CAMI1.K»*VU,K> 

CO  TO  172 
IC7  DO  I6S  IM.N 
DO  181  KM.N 

ICC  VNII  |.>N<I  |*EAN(1,RIIV(I,K) 

a  IF  A  IS  A  FUNCTION  OF  T,  ACaPORN  PCRtOOIC  OCCONAOStTlON.  t 
a  CcaUATE  AlTlDl.  and  ADJUST  PINAL  V(N*|)  a 

ititiiTiitttititttttiittaataaaattaaaaatatasstaaaaataiaiaaaataaaaaa 

ITD  IFdAT  .CO.  *)  CCTURN 
DO  12S  I«1.N 
DO  leS  J'l.N 
IK  C2F>«'  I.Jl'FIlD.I.JI 
17C  Th«T( I  )*N 

CALL  AFaTIAT.N.TNI 
DO  179  IM.N 
611  )>*.0* 

DO  179  JM.N 

173  611  )•dll*<AT(I.JI-A«I.JMaV(i.J) 

DO  174  IM.N 
0N|(4,4. I !•«.(* 


DO  174  J>t,N 

174  0Hi«4,4.i)>am(4.4,ii*ceAN(x,jia6(ji 
CO  17S  I'l.N 
DO  I7S  JM.N 

m  VN(Ii.VNIII-HtA2(I,J>tON|(4.4.J> 
at TURN 
M*  continue 


a  NonLiNrAa  nulti-2*$tcn  start*  ncrc 

a  ICOINNIfO  SECTION  DOCS  INITIALIZATION 


DO  24*  IM.N 

DO  24|  J«|.N 
DO  241  L«l.4 
PNI<L.J.II«*.M 
241  Pl(I.J)«*.D* 

VNIII**.0* 

IM  CONTInuC 

iPdNDCt  .or.  *)  00  TO  cc* 


a 

a 


mis  p.(i»  IS  SKI  auu.irr  "«MM»»>« 
PRO*  ^ 


l>AAOn 


I 

I 


TR 


t  tMO  MDC*  CM 
•  M  COtlTAlNS 


true  .or.  t>  M  TO  134 

:niT  .10.  Cl  ao  TO  zk 


I  CT  »Me  riMSM  0»  LOO"  2tS  — 

•  PMKl.I.Ji  CONTAINS  »NIC|,Oi  AMIII.I.J)  CONTAINS  ANKO.Ot 

•  UNAL  vin.ji  ARC  CALCULATCO  IN  LOC^  332 


DO  213  IM.N 
00  213  J«1.N 

I  Dtd.Jl— CAM(t.JI*UNXT(|.J> 

DO  21S  I'l.N 
DO  21S  J*1.N 
DO  21S  C«1.N 

'  DNI<3,|.Ji>AhI(3,|.J)-4N(I.K>«CAN(K.JI 
DO  311  1*1. N 

DO  3tD  J«t.M 

FCd.I.Jl'PHKI.t.JI 

DHI  ( 3, 1 .  J  »FN1 1 3, 1 .  J  >M1  d.  J  iM.OMaONd  .  J I 
)  Fe(2.I.JI*DHI(3.I.JI 
I  Dld.It’O.DO 
I  DO  2DC  I't.N 
DO  300  J«I.N 
FNICl.l.Jt'FCd.I.J) 

I  DM113.1. J>*FC(3,I,J> 

DO  330  KM.KSrcD 
CALI  CrmC.N.N.V.K.T.A) 

DO  331  1*1. N 
DO  331  J>I.N 

VNd  >*VNd  )*DNI(C.t.JU6(  JIOM 
I  CONTINUC 
DO  330  t'l.N 
DO  330  J«t.N 

I  D}d.l>>Pt«l.n-A3d.JltVN(J> 

DO  333  1*1. N 

IFdT  .CO.  0)  OSd.l>«O.M 
to  333  J'l.N 

I  DMKA.I.I  l•FHt<4,t.I)*CAM(X.3ltV(|.J) 

!  VMd»*Dld.l)*DMl«4,l.|) 

IFdAT  .CO.  0>  HCTUm 
DO  3S0  1*1. N 
to  390  J«t,N 
>  C3ANd.J»*^|(3,I.JI 
CO  TO  170 

0  3N9  Ottet  CM 

0  NCXT  BCLOU  307,  A3  CONTAINS  (ANIM(-3>  S 

I  eONTINu”********”***"*******"***********'”*****”***”**”**** 

IFdS  .CT.  tl  CO  TO  C7« 

IFdT  .CO.  01  CO  TO  330 

i  S?S  9^  FMXd.I^I  CONTAINS  FMXt3.L),L*0.1.C  • 

0  FINAL  V«N*|»  MC  CALCULATCD  IN  LOOF  333  00  83i  3 

K*270*X***«***************"******’*******”********************** 
to  370  J*l> 

Uld.X.J>«UNITlX,Jt-O.SOOtAN<X,J) 

yi<  3,  X ,  j  t'UNiTd .  j  )»o.sootAH<  X ,  j  > 

W3d.X.J)—3.B0tUNXrd.J>*AN3(i.J) 

M<3.X,J>*>3.30B(UNXr(|.J)«AHd.jM 

W3d,I.J>>Ullt.X.J> 

M3(|d,3>«UNXrtX.JI«l.S00SANd.J>«AN|<X.J> 

to  170  X*1#N 

DO  27?  J*t,N 

£0  274  L*1,N 

0Nl(|.t,Ji«OHXd.t.J)-Uld.l.L)«CAHlL.Jl 

aNI.3,I,J)«0Ml(B.I.J)-U3ll.I.liaEAH(l.Jl 

0Nl(3.t.J)«0Nl(3.I.J)-U3(l.l,Li«CAN«L.J> 

CONTINUC  ^  A 

QHI«l,I.J»>0MXt|,l.j>*wt(2,!.JI  .'vjp 

0HI«2.I.Jt*0Nl<t,l,J)*U3(f.X.Jt 
0H|(3.I.J>>0NI(3.X.Ji*«l3<t.X,JI 


374 
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c 

c 

c 

c 


e 

c 

c 

e 


rra  coHTiHuc 

<72  COMTINUE 

DO  27S  IM.3 
DO  27S  l<t,N 
DO  27S  JM.N 
DO  27S  K«1,N 

<7S  FICL.i.J).pHt(l,I,J) 

<7»  DO  2M  L*1.3 
DO  2D8  1*1. N 
DO  £D8  J>I.N 

<•< 

DO  KO  IM.H 
DO  2SO  K>1.3 
CAUL  ONtS.H.M.V.r.T.A) 

DO  ^9^  J-l.N 

vn(n<vN(ti*PMi(c.i,J>*c<jitH 
<t2  IFr(.E0.3>  VNI1>>VN>II«CAM(|,J»V(2.J) 
2M  C07<Tlri!JF 
CO  TO  293 
<tS  DO  2CS  1«I,N 
DO  23S  J>t.N 

<«S  VN!I>«VN>I)»EAH(I.J)BV(2,JI 
<93  IF<1AT  .EO.  •)  FCTuRtI 


•  FONCTIOH  of  T  PCfFOOn  FCRIOOIC  DCC0«9^S1T10M  • 

8  CVALU»TE  A(T(I>>.  and  adjust  FINAL  VIN»|»  • 

(ttnitltKlcsttBlfatmKBStltfKStMBtlFMBMMMDttnDDMttMM 

<97  TJ*T< 1  )«H 
T2»T»*H 

CALL  ArNT(C3EM,N,TI  I 
CALL  A^Mr<<.T,H,Ta> 

DO  294  I>|,N 
6<I >-9.D9 
fl(t.ll>9.C9 
DO  294  JM.N 

0(II>Ctn4«C3AH(I.JI'A(I.JM»V(<.Ji 
<94  Dt(t.l)>Ft(t.l)«(AT(l,J)>A(t.J>»Vt3,J) 

DO  23S  1«1,N 
Pt(2,II*«.D« 

DO  29S  J'l.N 

DO  2t«  1«|.N 

<9<  VNd  :«VNCII«HtPl«.|> 

HCTUtW 

<99  CO  2ii9  I't.N 
<99  VNtll>VN(I)>MtPl«,|> 
tCTU4N 
399  continue 
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C 

C 


c 

c 
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EAx'LxriHH),  AJ.iAMItll'J; 


«UP*K''EA 

IFilNftX  ,EO.  1 >  KUR*KSTCP*I 
DC  120  :>i.N 
329  VN<  I  CO. DR 

IF (IS  .CT.  tl  CO  TO  3<t 
IF( INDEX  .CO.  I >  CO  TO  357 
393  IF(IT  .EQ.  91  CO  TO  379 


•  CALCvLATC  PNI(3.K>.  C*9,I,2.  RESULTS  IN  AHKK.t.Jt  t 

00  3*4  1>!.N 
DO  394  J-1,N 

H1(1.I,J)*UNITII.J)*9.SIAN(I.J) 

MMa.I.Jl'UNlTU.Jl^l.SORtAHd.Jj^AHZlI.J) 

U2(I.I.J)— 2.C*«<UNIT<I.J)«AM(1,J)I 

U3<2,I.J>—2.De*UNIT(I.JI-4.D9tANCI,J)>3.09tAH2(I,J) 

M3<l.t.J)*WII2.I,JI 

394  lt3l<.l.J)*UNIT(l.J)*i.909tAN<I,JU3.M9AH2tl.J) 


mtt-wt  IS  BKi 

rRo*  oufY  pjs^vbhtr.  TO  jdc 
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M  3»*  l«l.N 

M  IM  J*I.N 

M  3*9  (•I.N 

.m.  41  M  to  SM 

0^4^J»0N|I  4.1.  J  IM(4(  1. 1.K  iOCONIK.  J I 

ir«KU4.C0.4l 

CdNTIMK 

CONVIWK 

8  310  KM.CUr 
310 

M  310  J*1.N 
00  310  L*I.O 

{;(t;OT».CO.O>  «t(l,i»OONllK.UJI 

iritNCCX.CO.OI  PKK.t.JfOMKC.l.}) 

IVilHDill.tO.tl  0MIt«.|,J>*OMl(K.Ljt-«N4«l.i»00Hl«KA.J> 

iriiNKM.co.ti  riiK.|,5i«OMi«K.|.}i 

M  311  >«I.KU0 
00  3tl  I«I.N 
to  3tl  JM.N 

IftlNOCK.CO.OI  0M|«.{.J»niK,I.J) 
iriiNrcM.co.t)  0Ni(0.i.ji*MiK.i.ji 
to  314  K*1.KU0 

£i*-i.y?‘***‘*''*'^***^*‘*» 

to  310  1*1,0 
00  310  3*1. N 

•atntattNMMnoimotnoMMtotttoottMooMooomuaaoowoooaoo 

a  MtcuuiTC  riiNki  v(N«ti  • 

aaaaatatatiaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 

VM(ii*v««<ti«NaON|«i,t,jiaf(j>  ^ 

iriK.co.Kuoi  «N(ii*Wi(ti4|aM<i.3>avc3.3» 

vOHTINUC 
00  TO  340 

00  m  1*1.0 

00  371  3*{.0 

70(1 1*^0111  >4<00<  1,3  »aV(3.3  • 
tOllOT  .n.  01  OgfuM 


ih:nd£.  .c j.  •<  GO  to  iot 

c  •4«tM««iitittitiiitta«fttii«*«(tt>fattatatt4titaitattiitaa44t44it 

C  t  T:  c|c  '5  •  a  f  jKflOt  or  T  • 

c  •••••t'tToatttMiaatitttintattataaaaitatataaaataaataaaaatatttttt 

Tl*T<  1  I4«< 

T3>1c»m 

COIL  OTNTlor.N.Tl I 
coil  OFNT«C1NM.O,TO| 

CM.I  aroTiri.o.Tli 
OO  341  t.i.N 
0(1 1*0.00 
Sit  (4. 3. 1 1*0.00 
001(4,4.11*0.00 
00  341  3*1.0 

«(ll*0(tt*(*t(t,J»^t|,j>)tV(l.3t 
001(4.3.1  |*0N|(4.3.|ltiC3llN(  1.3  >-0(1,31147(3.31 
341  0NI(4.4.t)*0NI(4.4.|l«iri(|.3>-O(t.3^iaV(4.3> 

00  340  t*|.N 
,  00  343  3*1.0 

•  • *  *0  f  0 

oCTuao 

c  oraaaaaaattaaaaaaaaataaataaaaaaiataaataaaattaaaaaataaiaaaaaaaiaata 

c  •  . .S**-®*-*^*  looiicjT  roi  ruNcrioo.  muj.ji,  3*o.i,t.3  a 

e  •••••■■■•••Mantraauaaaaaaaaataaaaaaaaaaiaaaaaaaiaiatoaaaaaoaaoa 

XT  iFdT  .10.  01  00  TO  370 

OO  3S0  1*1.0 
OO  3St  3*1.0 

01  ( 1 . 1 .3  K-ilOiri  { .  3 IFOMII  1,3 1/0.00 
01  «.|. 31* -uoir(  i.3)-aN(  1.1  >>ON0(i.3 1/3.00 
0i(l.t.3l*3.C0aUNlt(t.31««M(t.31-M(|.3l 


'J  V 
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M(«.  I.  J  1 .  J  (M.MaMHl  I.  J  I«1  t.  J ) 

|Mia^l^>UNlTI  I.)  )«a.NMIM  t.  J  )*11 . 


*UN  MMI 


i94CJ3SMt93f''99r-*S 
-9. 194$33s9saasi>««ac-9S 


9.CI  194'?46353424S9D»M 
•  .«I194?4I3144S*M34«9 


9. 1 39l9S94'»t499 1939*91 
9.139t9S9ST|S|6«999*91 


MX  9El  EMOR  •  9. 

M  •  9.63S99999C-9I 


9. 499191111 *99 


•9. 19493397749433910-99 
-9. 19493399999999990-99 


9.91 194759999973990*99 
9.91 19474931 4499990*99 


9. 13990993995174990*91 
9.13990999715199990*91 


MX  RCL  ERRIM  •  9.13310394994349980-99 

*  9.899999909*99  T  *  9.409990001*98 


-9. 19493390409798780-99 
-9 . 1 9453309989090990-95 


9.01104709909891970*99 

9.91194740314490090*99 


9.13009470930984710*91 

9.13000999719199990*91 


MX  RCL  ERROR  •  9.10599390437484840-99 

•  9.199099008*91  T  *  9.409999998*98 


•9. 19498a9488930S980-99 
-9.10493389909099900-95 


9.81199595340488700*99 
9.81 194748314489990*99 


9.13900989819793499*91 
9. 13000999715189990*91 


MX  REl  ERROR  •  R. 31486990891 398940-R4 

•  8.400000088*81  T  •  8.480808888*08 


•9. I844318908S076R30-08 
•9. 19493389800090900-99 


9.81 116789330988380*89 
9.81104748314489900*99 


9. 13008488814981780*01 
0.13000909719180088*91 


MX  RCL  ERROR  •  R. 58843991  i949188R0-93 
•  9.160009908*98  T  *  9.40i 


-9.10878818989354310-99 
•9. 19453380889100990-99 


TT4  — 


9.618W9I081391904D*99  9.13780158887950790*91 

9.81194748314489900*99  9.13889598715189990*91 


MX  881  ERROR  •  #.88741834483883589-98 
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